function [ results_fnames ] = osl_run_first_level_epoched_connectivity( oat )

% [ results_fnames ] = osl_run_first_level_epoched( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level trialwise GLM
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2011

global OSLDIR;

first_level=oat.first_level;
source_recon_name='source_recon';

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(~isfield(first_level,'bc_trialwise')),
    first_level.bc_trialwise=0;
end;

% check if user has asked to jailbreak the trialwise data
if(~isfield(first_level,'save_data'));
    first_level.save_data = 0;
end

if(~isfield(first_level, 'trialwise_directory'));
    first_level.trialwise_directory = [];
end

if(isfield(first_level, 'tf_time_downsample_factor'))
    multDSfactor  = 1./first_level.tf_time_downsample_factor; % user expresses time_downsample_factor as a divisive factor - as this is the default in SPM
else
    first_level.tf_time_downsample_factor = [];
end % if(~isempty(first_level.tf_time_downsample_factor)

if ~isempty(first_level.time_moving_av_win_size)
    warning('first_level.time_moving_av_win_size is not used with epoched data');
end;

contrast_list=first_level.contrast;
for c=1:length(contrast_list),
    contrast_list{c}=contrast_list{c}(:);
end;

clear stats;

if(sum(first_level.bc)>0 && first_level.baseline_timespan(1)>first_level.baseline_timespan(2)),
    error('first_level.baseline_timespan(1)>first_level.baseline_timespan(2)');
end;

if(~isfield(oat.first_level,'design_matrix') && ~isfield(oat.first_level,'design_matrix_summary')),
    error('Design matrix is not specified');
end;


for subi_todo=1:length(first_level.subjects_to_do),   
        
    sub=first_level.subjects_to_do(subi_todo);
      
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this subject
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D=source_recon_results.BF.data.D;
    if(~strcmp([oat.source_recon.dirname],D.path)),        
        D=path(D,[oat.source_recon.dirname]);
        D.save;
        D=spm_eeg_load([D.path '/' D.fname]);
        osl_save_oat_results(oat,source_recon_results);    
    end;
    
    %% setup things that are common to all subjects
    if(subi_todo==1),
              
      
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            error(['Connectivity not supported in sensor space']);
            
        else,
            
            % setup std space brain
            first_level_results.gridstep=source_recon_results.gridstep;
            
            S=[];
            S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
            S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
            S.current_level=first_level;
            S.lower_level_mni_coord=source_recon_results.mni_coord;
            S.lower_level_gridstep=source_recon_results.gridstep;
                    
            [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
            clear S;
        end;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up potential time frequency bases in tf_settings, note that any time or
        %% frequency averaging or sub-windowing will be done AFTER the
        %% time-frequency decomposition
        if(~strcmp(first_level.tf_method,'none')),
            isTFdata = 1;
            
            freq_range=source_recon_results.source_recon.freq_range; % range used for stat testing
            if(isempty(freq_range)),
                freq_range=first_level.freq_range;
                disp(['Using frequency range of ' num2str(freq_range) ' Hz']);
            end;        
      
            % setup Morlet Wavelet basis set
            if(strcmp(first_level.tf_method,'morlet')),
                disp('USING MORLET TF');
                % freq bins to use for morlet wavelets
                tf_frequencies=linspace(freq_range(1),freq_range(2),first_level.num_freqs);

                if(first_level.num_freqs==1)
                    tf_frequencies=(freq_range(1)+freq_range(2))/2;
                end;

                % setup wavelet basis set
                fres=D.fsample;
                tf_settings.morlet_basis = spm_eeg_morlet(first_level.tf_morlet_factor, 1000/fres, tf_frequencies);
                isTFdata = 1;

            elseif(strcmp(first_level.tf_method,'multitaper')),
                tf_frequencies=linspace(freq_range(1),freq_range(2),first_level.num_freqs);
                tf_settings.tf_freqres=freq_range(2)-freq_range(1);
                tf_settings.tf_timeres=0.05; %s

                first_level_results.tf_freqres=tf_settings.tf_freqres;
                isTFdata = 1;

            elseif(strcmp(first_level.tf_method,'hilbert')),
                disp('USING HILBERT TF');

                if(first_level.num_freqs>1)
                    tf_settings.tf_freqres=first_level.tf_hilbert_freq_res;
                    tf_frequencies=linspace(freq_range(1),freq_range(2),first_level.num_freqs);
                else
                    tf_frequencies=(freq_range(2)+freq_range(1))/2;
                    tf_settings.tf_freqres=freq_range(2)-freq_range(1);
                end;
                first_level_results.tf_freqres=tf_settings.tf_freqres;
                isTFdata = 1;
            else
                error('Invalid TF transform');
            end;
                
            if(isfield(tf_settings,'tf_freqres')),
                disp(['Freq bin width: ' num2str(tf_settings.tf_freqres)]);
            end;
            disp(['Freq bins are: ' num2str(tf_frequencies)]);       
                    
        else
            disp('IN TIME DOMAIN');
            % not doing time frequency
            tf_frequencies=-1;
            
        end;
        
        tf_settings.tf_method=first_level.tf_method;
        tf_settings.tf_logtransform=first_level.tf_logtransform;        
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
        
        if(isempty(first_level.time_range)),
           first_level.time_range=source_recon_results.source_recon.time_range;
        end;
        
        source_recon_time_indices=find(D.time>=source_recon_results.source_recon.time_range(1) & D.time<=source_recon_results.source_recon.time_range(2));
        source_recon_times=D.time(source_recon_time_indices);
        
        first_level_time_indices=find(source_recon_times>=first_level.time_range(1) & source_recon_times<=first_level.time_range(2));
        first_level_results.times=source_recon_times(first_level_time_indices);
        
        tres=first_level_results.times(2)-first_level_results.times(1);
        
        freq_indices=1;
        if(tf_frequencies>0)
            freq_indices=tf_frequencies>=first_level.freq_range(1) & tf_frequencies<=first_level.freq_range(2);
        end;
        first_level_results.frequencies=tf_frequencies(freq_indices);
        
        if(first_level.time_average),
            first_level_results.times=mean(first_level_results.times);
            tres=1;
        end;
        
        if(first_level.freq_average),
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end;
        
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);
        if(first_level.space_average),
            Nvoxels_out=1;
            first_level_results.mni_coord=mean(first_level_results.mni_coord,1);
        end;
            
        mni_coords=first_level_results.mni_coord;

        disp(['Using for seed: MNI coordinate ' num2str(first_level.connectivity_seed_mni_coord)]);

        % find the nearest index of the beamformed voxels to the specified mni_coord
        dists=(sqrt(sum((mni_coords-repmat(first_level.connectivity_seed_mni_coord,size(mni_coords,1),1)).^2,2)));

        [dist,first_level_results.seed_index]=min(dists);

        disp(['Nearest dipole for seed at MNI coordinate ' num2str(mni_coords(first_level_results.seed_index,:))]);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% do the following to set first_level_results.times                

        modality='MEG';
        chanind = strmatch(modality, D.chantype);
        chanind = setdiff(chanind, D.badchannels);
        if isempty(chanind)
            error(['No good ' modality ' channels were found.']);
        end
        
        weights=source_recon_results.BF.inverse.W.MEG{first_level_results.mask_indices_in_source_recon(first_level_results.seed_index)};
        
        if sum(isnan(weights))~=0,
            error('Seed has NaN weights');
        end;
                      
        dat_seed=zeros(1,length(source_recon_time_indices));
        dat_seed(1,:)=weights*D(chanind, source_recon_time_indices, 1);
        
        Sc=[];
        Sc.dat=dat_seed;
        Sc.dat_seed=dat_seed;
        Sc.times_in=source_recon_times;
        Sc.time_range=first_level.time_range;
        Sc.baseline_timespan=first_level.baseline_timespan;
        Sc.connectivity_seed_regress_zerolag=first_level.connectivity_seed_regress_zerolag;
        Sc.tf_time_moving_av_win_size=first_level.tf_time_moving_av_win_size;

        [conn, times_out]=feval(first_level.connectivity_method,Sc);
        first_level_results.times=times_out;
        first_level_time_indices=1:length(times_out);
        
        %%%%
                
        %%%%%%%%%%%%%%%%
        %% results containers
        ntpts=length(first_level_results.times);
        nfreqs=length(first_level_results.frequencies);        
        
        % do a quick check to determine the exact size of the dat array
        % after downsampling
        if ~isempty(first_level.time_downsample_factor)
            temp = [1:ntpts];
            temp = osl_spm_resample(temp,multDSfactor);

            ntptsGLM = size(temp,2);
            clear temp
        else
            ntptsGLM = ntpts;
        end
        first_level_results.stdcope=zeros([Nvoxels_out,ntptsGLM,length(contrast_list),nfreqs]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
        first_level_results.cope=zeros([Nvoxels_out,ntptsGLM,length(contrast_list),nfreqs]);% num_voxels x num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
        first_level_results.pseudo_zstat_var=zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)
    
    %%%%%%%%%%%%%%%%
    disp(['Reconstruct time courses and computing stats for dataset ' source_recon_results.fname]);   

    %% setup the GLM design matrix for this subject    
    if(isfield(first_level,'design_matrix')),
        first_level.design_matrix;
    else,    
        
        % setup trial types
        trialtypes=[];
        triallist=[];

        for i=1:length(source_recon_results.source_recon.conditions), % indexes conditions/triggers within subset
                               
            trigname=source_recon_results.source_recon.conditions{i};
            Ntrialspercond=length(D.pickconditions(trigname)); %% number of trials for this condition
            triallist=[triallist , D.pickconditions(trigname)];           
            trialtypes=[trialtypes ; ones(Ntrialspercond,1).*i];
            
        end;
        
        S3=[];
        S3.trialtypes=trialtypes;                
        
        if numel(first_level.design_matrix_summary)==length(oat.source_recon.results_fnames) &&(ischar(first_level.design_matrix_summary{sub})), 
          
            S3.Xsummary=first_level.design_matrix_summary{sub};
        
            if(isfield(first_level,'trial_rejects'))
                S3.trial_rejects=first_level.trial_rejects;
            else
                error('Need to set oat.first_level.trial_rejects');
            end;
        else
            S3.Xsummary=first_level.design_matrix_summary;
        end;        
 
        first_level.x=setup_beamformer_designmatrix(S3);
        
        
        const_reg=0;        
        for p=1:size(first_level.x,2),
            if(range(first_level.x(:,p))==0),
                if(const_reg),
                    warning('Multiple constant regressors in the design matrix!!!!?');
                end;
                const_reg=1;
                
            end;                                    
        end;
        
        if(const_reg && first_level.do_glm_demean)
            warning('CAREFUL: are you sure you want the first_level.do_glm_demean flag on with constant regressors in the design matrix!!!!?');
        elseif(~const_reg && ~first_level.do_glm_demean)
            warning('CAREFUL: are you sure you want no first_level.do_glm_demean flag on with no constant regressors in the design matrix!!!!?');
        end;       
        
    end;
    
    Ntrials   = length(triallist);

    if ~isempty(first_level.tf_time_downsample_factor);
    	first_level.x = osl_spm_resample(first_level.x,multDSfactor);
    end;
    
    % precompute some things
    contrast_list=first_level.contrast;
    x=first_level.x;
    pinvxtx=pinv(x'*x);
    pinvx=pinv(x);
    
    %%%%%%%%%%%%%%%
    %% setup channels
    modality='MEG';
    chanindmeg = strmatch(modality, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=ones(length(chanindmeg),1);
        chanind = chanindmeg;
    end; 
        
    %%%%%%%%%%%%%%%
    %% load in sensor data 
    first_level_results.D_sensor_data=D;
    sensor_data=D(chanind, source_recon_time_indices, triallist);
    first_level_results.chanind=chanind;

    %%%%%%%%%%%%%%%%
    %% connectivity
    % recon seed time course        
    weights=source_recon_results.BF.inverse.W.MEG{first_level_results.mask_indices_in_source_recon(first_level_results.seed_index)};
        
    if sum(isnan(weights))~=0,
        error('Seed has NaN weights');
    end;

    dat_seed=zeros(Ntrials,length(source_recon_time_indices));
    for tri=1:Ntrials, % indexes trials
        dat_seed(tri,:)=weights*sensor_data(:, :, tri);
    end;
            
    if ~isempty(first_level.tf_time_downsample_factor);
        disp(['Downsampling by a factor of ' num2str(first_level.tf_time_downsample_factor) ' for the GLM (data at ' num2str(1/tres) 'Hz; GLM at ' num2str((1/tres)./first_level.tf_time_downsample_factor) 'Hz).']);
    else
        disp('Not downsampling for the GLM');
    end
    disp('First level COPEs outputted will have dimension Nvoxels x Ntpts x Ncontrasts x Nfreqs:');
    tmp=[size(first_level_results.cope) 1 1];
    disp([num2str(tmp(1:4))]);
        
    %%%%%%%%%%%%%%%%
    %% do everything from now on on one voxel at a time (for the sake of
    %% RAM) in the first level mask space    
    if first_level.save_data
        % do a quick check to determine the exact size of the dat array
        % after downsampling
        if ~isempty(first_level.tf_time_downsample_factor);
            temp = [1:Nsamples];
            temp = osl_spm_resample(temp,multDSfactor);
            NsamplesTrlwise = size(temp,2);
            clear temp
        else
            NsamplesTrlwise = Nsamples;
        end
        
        if isTFdata
            trlwise.type    = first_level.tf_method;
            trlwise.dat     = nan(Nvoxels_multidipole,size(source_recon_results.Y_evoked,1),Ntrials,NsamplesTrlwise,nfreqs); % if not a tf, the last dimenstion will be 'vestigial'
            trlwise.freqs   = first_level_results.frequencies;
            trlwise.times   = linspace(first_level_results.times(1), first_level_results.times(end), NsamplesTrlwise);
            
        else
            trlwise.type = 'time_domain';
            % trlwise dat will be assigned below
        end
        trlwise.mask        = read_avw([oat.source_recon.dirname '/' oat.first_level.name '_mask']);
        trlwise.mnicoords   = first_level_results.mni_coord;
        trlwise.gridstep    = first_level_results.gridstep;
    end
    
    %%%%%%%%%%%%%%%%
    %% now do space averaging 
    if(first_level.space_average),
        
        dat=zeros(Ntrials,length(source_recon_time_indices));
        num=0;
        
        for indind=1:length(first_level_results.mask_indices_in_source_recon), % indexes brain space
            if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis                          
                weights=zeros(1,length(chanind));
                % for sensor space - just compute for all sensors
                weights(indind)=1;
            else
                weights=source_recon_results.BF.inverse.W.MEG{first_level_results.mask_indices_in_source_recon(indind)};            
                first_level_results.pseudo_zstat_var(indind)=weights*weights';
                weights=weights/sqrt(first_level_results.pseudo_zstat_var(indind)); % to compute pseudo z-stats
            
            end;
            
            if sum(isnan(squash(weights)))==0,    
                for tri=1:Ntrials, % indexes trials
                    dat(tri,:)=dat(tri,:)+weights*sensor_data(:, :, tri);
                end;
                num=num+1;
            end;
        end;
        dat=dat/num;
    end;

    tfdat_seed=[];
    ft_progress('init', 'etf');

    for indind=1:Nvoxels_out, % indexes brain space
            
        ft_progress(indind/Nvoxels_out);                                   

        if(~first_level.space_average),
            if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis                          
                weights=zeros(1,length(chanind));
                % for sensor space - just compute for all sensors
                weights(indind)=1;
            else
                weights=source_recon_results.BF.inverse.W.MEG{first_level_results.mask_indices_in_source_recon(indind)};            
                first_level_results.pseudo_zstat_var(indind)=weights*weights';
                weights=weights/sqrt(first_level_results.pseudo_zstat_var(indind)); % to compute pseudo z-stats
             end;

            if sum(isnan(squash(weights)))==0, 
                dat=zeros(Ntrials,length(source_recon_time_indices));
                for tri=1:Ntrials, % indexes trials
                    dat(tri,:)=weights*sensor_data(:, :, tri);
                end;
            else
                dat=nan;
            end;
        end;

        if sum(isnan(squash(dat)))~=0,

            first_level_results.cope(indind,:,:,:)=nan;
            first_level_results.stdcope(indind,:,:,:)=nan;

        else,        

            Sc=[];
            Sc.dat=dat;
            Sc.dat_seed=dat_seed;
            Sc.times_in=source_recon_times;
            Sc.time_range=first_level.time_range;
            Sc.baseline_timespan=first_level.baseline_timespan;
            if ~isempty(tfdat_seed),
                Sc.tfdat_seed=tfdat_seed;
            end;
            Sc.tf_time_moving_av_win_size=first_level.tf_time_moving_av_win_size;            
            Sc.connectivity_seed_regress_zerolag=first_level.connectivity_seed_regress_zerolag;

            [conn, first_level_results.times, tfdat_seed]=feval(first_level.connectivity_method,Sc);
            dat=conn;

            if(first_level.bc_trialwise),
                baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);

                if(sum(baseline_time_indices)>0)
                    % normalise by noise in baseline window:
                    for tri=1:size(dat,1), % indexes trials
                        for f = 1 : size(dat,3),
                            dat(tri,:,f)=(dat(tri,:,f)-mean(dat(tri,baseline_time_indices,f)));
                        end;
                    end;
                else
                    warning('Could not baseline correct as there are no timepoints in the specified baseline period');
                end;
            end;

            %% now do windowing over time and frequency
            dat=dat(:,first_level_time_indices,freq_indices);

            %% now do averaging over time and frequency
            if(first_level.time_average)
                dat=mean(dat,2);
            end;

            if(first_level.freq_average)
                dat=mean(dat,3);
            end;

            %% do temporal downsampling if required

            % downsample in time
            if ~isempty(first_level.tf_time_downsample_factor),

                dat = osl_spm_resample(dat,multDSfactor);
                % correct the timebase
                first_level_results.times = linspace(min(first_level_results.times),max(first_level_results.times),size(dat,2));
                % resize the results containers
                first_level_results.cope    = first_level_results.cope(:,1:size(dat,2),:,:,:);
                first_level_results.stdcope = first_level_results.stdcope(:,1:size(dat,2),:,:,:);

                if isTFdata
                    % check for <=0s in new dat.  Give a warning if there are any
                    if min(min(dat)) <= 0
                        warning('There are zeros or negative values in the downsampled time-frequency data! The stdcope will be complex...');
                    end % if min(min(dat)) <= 0
                end % if isTFdata

            end % if ~isempty(first_level.tf_time_downsample_factor);


            %% Do GLM
            if first_level.doGLM
                if(size(dat,1)>1)
                    %% now do the GLM stats
                    for f = 1 : size(dat,3),
                        for t= 1 : size(dat,2), % indexes time points within trials

                            if(first_level.use_robust_glm),
                                [b,stats2]=robustfit(x,(dat(:,t,f)),'bisquare',4.685,'off'); % fits GLM robustly, Y=XB+E
                                for c=1:length(contrast_list),

                                    if(strcmp(oat.first_level.cope_type,'coape')),
                                        first_level_results.cope(indind,t,c,f)=contrast_list{c}'*abs(b);
                                    elseif(strcmp(oat.first_level.cope_type,'cope')),
                                        first_level_results.cope(indind,t,c,f)=contrast_list{c}'*b;
                                    elseif(strcmp(oat.first_level.cope_type,'acope')),
                                        first_level_results.cope(indind,t,c,f)=abs(contrast_list{c}'*b);
                                    else
                                        error('Invalid cope_type');
                                    end;
                                    first_level_results.stdcope(indind,t,c,f)=sqrt(contrast_list{c}'*stats2.covb*contrast_list{c});
                                end;
                                pe(t,:)=b; % num_timepoints x  num_pes
                            else
                                [copeout, varcopeout, coapeout, dof,peout]=glm_fast_for_meg((dat(:,t,f)),x,pinvxtx,pinvx,contrast_list,0);
                                if(strcmp(oat.first_level.cope_type,'coape')),
                                    first_level_results.cope(indind,t,:,f)=coapeout;
                                elseif(strcmp(oat.first_level.cope_type,'cope')),
                                    first_level_results.cope(indind,t,:,f)=copeout;
                                elseif(strcmp(oat.first_level.cope_type,'acope')),
                                    first_level_results.cope(indind,t,:,f)=abs(copeout);
                                else
                                    error('Invalid cope_type');
                                end;

                                first_level_results.stdcope(indind,t,:,f)=sqrt(varcopeout);
                                pe(t,:)=peout; % num_timepoints x  num_pes
                            end;

                        end;

                        %% do bc on the copes
                        if(sum(first_level.bc)>0),
                            baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);

                            for c=1:length(contrast_list),
                                if(first_level.bc(c))
                                    first_level_results.cope(indind,:,c,f)=first_level_results.cope(indind,:,c,f)-mean(first_level_results.cope(indind,baseline_time_indices,c,f),2);
                                end;
                            end;

                        end;
                    end;
                end;
            else
                first_level_results.cope(indind,:,1,:)=dat(1,:,:);
                first_level_results.stdcope(indind,:,1,:)=1;
            end;
        end;

        first_level_results.cope(first_level_results.seed_index,:,:,:)=0;
        first_level_results.stdcope(first_level_results.seed_index,:,:,:)=1e16;    

    end;

    first_level_results.cope(first_level_results.seed_index,:,:,:)=1;
    first_level_results.stdcope(first_level_results.seed_index,:,:,:)=1e16;

    ft_progress('close');
    
    %% save results
    first_level_results.x=x;
    first_level_results.contrasts=first_level.contrast;
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.S=first_level;
    first_level_results.bf_S=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    first_level_results.subject_name=source_recon_results.subject_name;
    
    % save results to disk
    if first_level.save_data
        if ~isempty(first_level.trialwise_directory)
            trlWise_fname=fullfile(first_level.rawTF_directory, [first_level_results.subject_name '_trlWise']);
        else
            first_level_results.fname=[ first_level_results.subject_name '_' first_level_results.name ];
            trlWise_fname=fullfile(oat.source_recon.dirname, [first_level_results.subject_name '_trlWise']);
        end
        disp(['Saving trialwise data in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        save(trlWise_fname, 'trlwise', '-v7.3');
        clear trlwise
    end % if first_level.save_data
    
    if first_level.doGLM
        first_level_results.fname=[ first_level_results.subject_name '_' first_level_results.name ];
        disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        osl_save_oat_results( oat, first_level_results );
        results_fnames{sub}=first_level_results.fname;
    end % if first_level.doGLM
    
    if (strcmp(source_recon_results.recon_method,'none') | first_level.space_average),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
end;


